clear all
set more off
cd /Users/zimaoxiao/Dropbox/CC_Yield_Predict/replication_package/  // Navigate to replication folder on your own machine
global weather "gdd10_29 gdd29 prec prec2"

qui {
* ============================
* 	  Linear trends spec.
* ============================
use if inrange(year,1950,2000) using data/regression/reg_2020_prediction_gdd, clear
	
* apply trends and yield data filter
merge m:1 county_fips using data/yield/filter, keep(3) nogenerate
	
* regression
reghdfe lncornyield $weather state_fips#c.trend [aweight=corn_areaHarv], absorb(county_fips) vce(cluster county_fips StateXYear)

* save the estimation results
tempfile corn_st_lin
parmest, label saving(`"`corn_st_lin'"',replace) format(p %8.2f) stars(0.1 0.05 0.01)

* extract temperature coefficients		
use parm estimate using `corn_st_lin', clear
keep if _n<=2
foreach i in gdd10_29 gdd29 {
	gen `i' = estimate if parm == "`i'"
	summ `i'
	local k = r(mean)
	replace `i' = `k'
}
keep if _n==1
drop parm estimate
tempfile corntempcoef
save `corntempcoef', replace

* extract county linear trends coefficients		
use parm estimate using `corn_st_lin', clear // change "30more_unavailable" to change filter	
drop if _n<5
gen quad = ustrright(parm,1) // linear trend = "d", quad trend = "2"
drop if quad != "d"
drop quad
destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
drop parm			
rename estimate stlin
tempfile cornlintrend
save `cornlintrend', replace
						
* put coefficients and calculated climate change together
use county_fips CC* using data/regression/reg_2020_prediction_gdd, clear
bys county_fips: keep if _n == 1
gen state_fips = floor(county_fips/1000)
merge m:1 state_fips using `cornlintrend', keep(3) noreport nogenerate
cross using `corntempcoef'
	
* calculate climate change impacts by county
bys county_fips: gen ccimpact_corn = (CC_gdd10_29*gdd10_29) + (CC_gdd29*gdd29)
								
* calculate overall impacts by county (climate change + tech effects)
* 1950 as the 1st year, t=1; 2000 as the 51st year, t=51; 2020 as the 71st year, t=71
bys county_fips: gen trends_corn = (CC_gdd10_29*gdd10_29) + (CC_gdd29*gdd29) + (stlin*20)
keep county_fips ccimpact_corn trends_corn 
tempfile corn_cccoef
save `corn_cccoef', replace
	
* merge in regression coefficients, calculate predicted 2020 yields
use data/yield/filter, clear
merge 1:1 county_fips using `corn_cccoef', keep(3) noreport nogenerate
	
* predicted 2020 yields: climate change
bys county_fips: gen cc2020_lin = cornyield_2000*(1+ccimpact_corn)

* predicted 2020 yields: climate change + tech effects
bys county_fips: gen trends2020_lin = cornyield_2000*(1+trends_corn)

* Table S2: rows (1), (2), (4) and (6)		
rename (cornyield_2000 cornyield_2020) (yield2000 yield2020)
keep county_fips yield2020 yield2000 cc2020_lin trends2020_lin
tempfile gdd_2020_lin
save `gdd_2020_lin', replace

* ============================
* 	 Quadratic trends spec.
* ============================
use if inrange(year,1950,2000) using data/regression/reg_2020_prediction_gdd, clear
	
* apply trends and yield data filter
merge m:1 county_fips using data/yield/filter, keep(3) nogenerate
				
* regression
reghdfe lncornyield $weather state_fips#c.trend state_fips#c.trend2 [aweight=corn_areaHarv], absorb(county_fips) vce(cluster county_fips StateXYear)

* save the estimation results
tempfile corn_st_quad
parmest, label saving(`"`corn_st_quad'"',replace) format(p %8.2f) stars(0.1 0.05 0.01)
											
* extract temperature coefficients		
use parm estimate using `corn_st_quad', clear
keep if _n<=2
foreach i in gdd10_29 gdd29 {
	gen `i' = estimate if parm == "`i'"
	summ `i'
	local k = r(mean)
	replace `i' = `k'
}
keep if _n==1
drop parm estimate
tempfile corntempcoef
save `corntempcoef', replace
						
* extract county linear trends coefficients		
use parm estimate using `corn_st_quad', clear
drop if _n<5
gen quad = ustrright(parm,1)
preserve
	* first-order
	drop if quad != "d"
	drop quad
	destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
	drop parm			
	rename estimate stlin
	tempfile cornlintrend
	save `cornlintrend', replace
restore
	* second-order
	keep if quad == "2"
	destring parm, generate(state_fips) ignore("b.state_fips#c.trend")
	gen state_fips_1 = floor(state_fips/10)
	replace state_fips = state_fips_1
	drop state_fips_1
	drop parm quad
	rename estimate stquad
	tempfile cornquadtrend
	save `cornquadtrend', replace
	
* put coefficients and calculated climate change together
use county_fips CC* using data/regression/reg_2020_prediction_gdd, clear
bys county_fips: keep if _n==1
gen state_fips = floor(county_fips/1000)
merge m:1 state_fips using `cornlintrend', keep(3) noreport nogenerate
merge m:1 state_fips using `cornquadtrend', assert(3) keep(3) noreport nogenerate
cross using `corntempcoef'
	
* calculate climate change impacts by county
bys county_fips: gen ccimpact_corn = (CC_gdd10_29*gdd10_29) + (CC_gdd29*gdd29)
				
* calculate overall impacts by county (climate change + tech effects)
* 1950 as the 1st year, t=1; 2000 as the 51st year, t=51; 2020 as the 71st year, t=71
bys county_fips: gen trends_corn = (CC_gdd10_29*gdd10_29) + (CC_gdd29*gdd29) + (stlin*20) + (stquad*2440)
keep county_fips ccimpact_corn trends_corn 
tempfile corn_cccoef
save `corn_cccoef', replace
	
* merge in regression coefficients, calculate predicted 2020 yields
use data/yield/filter, clear
merge 1:1 county_fips using `corn_cccoef', keep(3) noreport nogenerate
				
* predicted 2020 yields: climate change
bys county_fips: gen cc2020_quad = cornyield_2000*(1+ccimpact_corn)
				
* predicted 2020 yields: climate change + tech effects
bys county_fips: gen trends2020_quad = cornyield_2000*(1+trends_corn)

* Table S2: rows (3) and (5)		
rename (cornyield_2000 cornyield_2020) (yield2000 yield2020)
keep county_fips cc2020_quad trends2020_quad
tempfile gdd_2020_quad
save `gdd_2020_quad', replace

* ============================
* 	         Output
* ============================
use `gdd_2020_lin', clear
merge 1:1 county_fips using `gdd_2020_quad', assert(3) keep(3) nogenerate
noi logout, save(output/TableS5) word dec(1) replace: tabstat yield2000 cc2020_lin cc2020_quad trends2020_lin trends2020_quad yield2020, stat(p10 p25 p50 p75 p90) format(%10.4f) c(s)
save data/figure/gdd_2020, replace // for Figure 2

}
*** EOF
